#' @TODO 生成pdf、-72.tif、-300.tif三张图
#' @title 给定ggplot绘图对象p，生成pdf、-72.tif、-300.tif三张图
#' @param od 图片输出路径
#' @param name 图片名称
#' @param w 图片宽度，默认6
#' @param h 图片高度，默认6
#' @param p 图片存储对象
#' @export 
#' @return NULL
#' @examples plotout(p = p, num = 1, od = 'out_home')
#' @author *CY*
#'
plotout <- function(
    p = NULL, od = NULL, name = NULL, w = 6, h = 6, num = NULL,
    savetiff = F, plot_tif = savetiff) {
  if (!is.null(num)) {
    name <- num
  }
  if (!dir.exists(od)) {
    dir.create(od, recursive = TRUE)
  }
  pdf(
    file = paste0(od, "/Figure_", name, ".pdf"), width = w,
    height = h
  )
  print(p)
  dev.off()
  if (isTRUE(savetiff)) {
    tiff(
      file = paste0(od, "/Figure_", name, "-72.tif"),
      width = w, height = h, units = "in", res = 72
    )
    print(p)
    dev.off()
    tiff(
      file = paste0(od, "/Figure_", name, "-300.tif"),
      width = w, height = h, units = "in", res = 300
    )
    print(p)
    dev.off()
    png(
      file = paste0(od, "/Figure_", name, "-300.png"),
      width = w, height = h, units = "in", res = 300
    )
    print(p)
    dev.off()
  }
}

#'
#' @TODO ## 差异分析
#' @title 差异分析
#' @description 利用R包`limma`进行差异分析。输入表达谱与分组数据，进行差异分析。使用标准化后的数据进行差异分析。
#' @param exp 表达谱数据
#' @param group_infor 分组信息, 分组信息列名以`group`做列名,样本所在列列名为`sample`
#'  > 注： group中默认第一会作为`control`，另一分组作为`case`。
#' 例如 tumor与normal，字母n在t前面，所以第一顺位是normal，以normal 为对照，观察tumor的上调与下调
#' @param output_dir 输出文件目录
#' @param var_name 用于文件命名
#' @param levels_manual 指定差异分析中的levels，如果是NULL，则默认group_infor中第一顺位的分组为对照组，第二个为实验组。
#' 否则指定一个顺序，是group_infor中的组别名称。eg：levels_manual = c('Normal','Tumor')
#'
#' @return  a *list*，第一个元素全部degs结果，第二个为表达谱，是按分组文件筛选后的表达谱
#' @export
#'
#' @author *WYK*
#'
DEGs <- function(exp = NULL, group_infor = NULL, output_dir = "./", var_name = NULL, levels_manual = NULL) {
  if (is.null(levels_manual)) {
    levels_manual <- c(sort(unique(group_infor$group))[1], sort(unique(group_infor$group))[2])
  }

  if (length(var_name) == 0) {
    var_name <- paste0(sample(c(letters, LETTERS), 5), collapse = "")
  }

  suppressMessages(library(magrittr))

  library(tidyverse)
  group_infor <- na.omit(group_infor)
  rownames(group_infor) <- NULL

  if (!"group" %>% is_in(colnames(group_infor))) {
    warning("缺少分组信息列，或者 分组信息所在列列名不是group，程序跳出。")
    return()
  }

  common_sample <- intersect(group_infor$sample, colnames(exp))
  group_infor <- group_infor %>%
    column_to_rownames("sample") %>%
    .[common_sample, , drop = F] %>%
    rownames_to_column("sample")
  exp_1 <- exp[, common_sample]

  colnames(exp_1) == group_infor$sample

  library(limma)
  eset <- as.matrix(exp_1)

  group <- factor(group_infor$group, levels = levels_manual, ordered = F)
  design <- model.matrix(~group)
  colnames(design) <- levels(group)

  fit <- lmFit(eset, design)
  fit2 <- eBayes(fit)
  dif <- topTable(fit2, coef = 2, n = nrow(fit2))
  # head(dif)
  degs <- dif %>%
    as.data.frame() %>%
    rownames_to_column(var = "gene") %>%
    dplyr::rename(log2FC = logFC, pvalue = P.Value, adjusted_pvalue = adj.P.Val) %>%
    dplyr::select(gene, log2FC, pvalue, adjusted_pvalue)

  if (!dir.exists(sprintf("%s/", output_dir))) {
    dir.create(sprintf("%s/", output_dir), recursive = T)
  }

  degs %>% write_tsv(
    x = ., file = sprintf("%s/deg_res_%s.txt", output_dir, var_name),
    quote = "none",
  )

  return(degs)
}

#' @title  ## 按阈值筛选差异分析结果
#' @TODO  按阈值筛选差异分析结果
#' @description 按设定阈值对差异分析结果进行分析。差异分析全部结果与各项阈值.
#' @param degs gene、log2FC等数据
#' @param log2FC_value log2FC绝对值，默认1，即两倍差异
#' @param DEG_adjusted_pvalue 矫正后p值筛选结果，默认0.05
#' @param DEG_pvalue 如果`DEG_adjusted_pvalue`不支持后续分析，可以使用p值筛选结果。
#'   > 注：如果`DEG_pvalue != NULL`，在程序运行中，其优先级高于`DEG_adjusted_pvalue`，
#'   >     即`DEG_pvalue`与`DEG_adjusted_pvalue`都存在的情况下，会优先使用`DEG_pvalue`
#' @param intrest_gene 兴趣基因集，字符串向量,如果为NULL，则不对差异结果进行筛选
#'
#' @export
#' @author *WYK*
#'
model_data <- function(degs = NULL, log2FC_value = 1, DEG_adjusted_pvalue = 0.05,
                       DEG_pvalue = NULL, intrest_gene = NULL) {
  if (is.null(DEG_pvalue)) {
    deg1 <- degs %>%
      filter(adjusted_pvalue < DEG_adjusted_pvalue) %>%
      filter(log2FC > log2FC_value | log2FC < (-log2FC_value))
  } else {
    deg1 <- degs %>%
      filter(pvalue < DEG_pvalue) %>%
      filter(log2FC > log2FC_value | log2FC < (-log2FC_value))
  }

  up <- deg1 %>%
    filter(log2FC > log2FC_value) %>%
    nrow(.)
  down <- deg1 %>%
    filter(log2FC < -log2FC_value) %>%
    nrow(.)


  if (is.null(DEG_pvalue)) {
    msg_tmp <- sprintf(
      "There're %s  DEgenes with log2FC = %s, adj.p = %s. Up:%s , Down:%s .", nrow(deg1),
      log2FC_value, DEG_adjusted_pvalue, up, down
    )
    print(msg_tmp)
  } else {
    msg_tmp <- sprintf(
      "There're %s  DEgenes with log2FC = %s, p = %s. Up:%s , Down:%s .", nrow(deg1),
      log2FC_value, DEG_pvalue, up, down
    )
    print(msg_tmp)
  }


  if (length(intrest_gene != 0)) {
    deg2 <- deg1[which(deg1$gene %in% intrest_gene), ]
    if (nrow(deg2) == 0) {
      sprintf(
        "In \"%s\"'s setting, 0 gene mapped in gene list.",
        msg_tmp
      )
    } else {
      msg_tmp_1 <- sprintf("In \"%s\"'s setting, %s gene mapped in gene list.", msg_tmp, nrow(deg2))
      message(msg_tmp_1)
    }
  } else {
    deg2 <- deg1
  }

  # tmp <- list(deg2 = deg2, exp_2 = exp_2)

  return(deg2)
}

#' @title  将p值转化为星号
#' @description 将p值转化为星号，结果
#' @param p 数字型向量
#' @return 返回一个data.frame，第一列为p值，第二列为星号
#' @author *WYK*
#'
formatp <- function(p) {
  stopifnot(!class(p) %in% c("character", "list"))
  p_sig <- cut(p, c(0, .01, .05, Inf), c("**", "*", ""), include.lowest = TRUE)
  dplyr::tibble(p, p_sig)
}

#' @title  ## 绘制差异基因表达热图
#' @description 按照degs结果来绘制热图，自动拟合图形宽高
#' @param exprs 表达谱
#' @param group 分组信息dataframe，两列，第一列sample样本列，第二列group分组新系列
#' @param degs 差异分析结果
#' @param od 结果输出路径
#' @return NULL
#' @export
#' @author *WYK*
#'
deg_heat <- function(
    exprs = dat$data_exprs, group = dat$group,
    degs = xgne_in_deg, p_type = "pvalue", od = NULL,annotation_colors = NULL) {
  
  unique(group[[2]])

  degs$p_sig <- formatp(degs[[p_type]])[[2]]
  degs$symbol <- str_c(degs$gene, " ", degs$p_sig)

  df_for_heat <- exprs[degs$gene, ] %>%
    rownames_to_column("gene") %>%
    pivot_longer(-gene, values_to = "exprs", names_to = "sample") %>%
    inner_join(group) %>%
    inner_join(degs %>% select(gene, symbol)) %>%
    arrange(group, exprs)

  p <- tidyheatmaps::tidy_heatmap(df_for_heat,
    rows = "symbol", values = "exprs", columns = "sample",
    annotation_col = "group", scale = "row",annotation_colors = annotation_colors,
    show_colnames = F, cluster_rows = T, color_legend_min = -2, color_legend_max = 2,
    colors = c("#0C6291", "white", "#A63446"), silent = T
  )

  d <- data.frame(x = c(14, 10, 6, 4), y = c(4.75, 3.5, 3.25, 3.25))
  fit <- lm(y ~ x, data = d)
  p_length <- predict(fit, tibble(x = length(degs[[1]])))

  plotout(p = p, name = "deg_heat", od = od, w = 5.2, h = p_length)
}


#' @title 多变量同一分组在一个boxplot中
#' @description 多变量同一分组在一个boxplot中
#' @param input data.frame，输入数据
#' @param x_ 字符串向量，x轴坐标对应列
#' @param y_ 字符串向量，y轴坐标对应列
#' @param group_ 字符串向量，分组信息对应列
#' @param type_ y轴数值类型，比如说
#' @param color_in_p 字符串向量，分组对应颜色
#' @param x_angle 数字，x轴字体倾斜角度
#' @param style 单一字符串，可选box1，box2，vln。box1：白线中线箱线图；box2：填充半透明箱线图；vln：填充半透明的小提琴图
#' @export
#' @return ggplot对象
#' @author *WYK*
#'
Box_in_One_p <- function(
    input, x_, y_, group_, color_in_p = RColorBrewer::brewer.pal(4, "Set1"),
    x_angle = 45, style = "box1") {
  library(tidyverse)

  ignoreCase <- c(
    "NA", "NAN", "NaN", "NX", "MX", "TX", "Not Reported", "GX", "UNK", "-", "", " ",
    "unknown", "UNknown", "notreported", "gx", NA
  )

  input <- input %>%
    dplyr::filter(!group_ %in% ignoreCase)

  color_in_p <- color_in_p[seq_len(length(unique(input[[group_]])))]

  p4 <- switch(style,
    box1 = {
      ggplot(
        data = input,
        aes_string(x_, y_, fill = group_)
      ) +
        geom_boxplot(aes_string(col = group_),
          outlier.size = .2, outlier.shape = NA,
          width = .5,
          position = position_dodge(width = .63)
        ) +
        xlab(NULL) +
        ggpubr::theme_pubr(13) +
        theme(
          axis.text.x = element_text(angle = x_angle, hjust = 1, vjust = 1, size = 10),
          legend.position = "bottom",
          legend.title = element_blank()
        ) +
        ggpubr::stat_compare_means(
          aes(
            group = get0(group_),
            label = after_stat(p.signif)
          )
        ) +
        ggplot2::coord_cartesian(clip = "off") +
        scale_fill_manual(values = color_in_p) +
        scale_color_manual(values = color_in_p) +
        stat_summary(
          fun = median, geom = "crossbar", color = "white", width = 0.59, size = .2,
          position = position_dodge(width = .63)
        )
    },
    box2 = {
      ggplot(
        data = input,
        aes_string(x_, y_, fill = group_)
      ) +
        geom_boxplot(aes_string(col = group_),
          outlier.size = .2, outlier.shape = NA,
          width = .5,
          position = position_dodge(width = .63)
        ) +
        xlab(NULL) +
        ggpubr::theme_pubr(13) +
        theme(
          axis.text.x = element_text(angle = x_angle, hjust = 1, vjust = 1, size = 10),
          legend.position = "bottom",
          legend.title = element_blank()
        ) +
        ggpubr::stat_compare_means(
          aes(
            group = get0(group_),
            label = after_stat(p.signif)
          )
        ) +
        ggplot2::coord_cartesian(clip = "off") +
        scale_fill_manual(values = ggplot2::alpha(color_in_p, .75)) +
        scale_color_manual(values = color_in_p)
    },
    vln = {
      ggplot(
        data = input,
        aes_string(x_, y_, fill = group_)
      ) +
        # ggplot2::geom_violin()
        geom_violin(aes_string(col = group_),
          width = .5,
        ) +
        scale_fill_manual(values = ggplot2::alpha(color_in_p, .75)) +
        scale_color_manual(values = color_in_p) +
        xlab(NULL) +
        ggpubr::theme_pubr(13) +
        theme(
          axis.text.x = element_text(angle = x_angle, hjust = 1, vjust = 1, size = 10),
          legend.position = "bottom",
          legend.title = element_blank()
        ) +
        ggpubr::stat_compare_means(aes(group = get0(group_), label = after_stat(p.signif))) +
        ggplot2::coord_cartesian(clip = "off")
    },
    mix1 = {
      ggplot(
        data = input,
        aes_string(x_, y_, fill = group_)
      ) +
        geom_violin(alpha = .8, width = .7) +
        geom_boxplot(fill = "white", width = .1, outlier.shape = NA) +
        xlab(NULL) +
        ggpubr::theme_pubr() +
        scale_fill_manual(values = ggplot2::alpha(color_in_p, .75)) +
        ggpubr::stat_compare_means(aes(group = get0(group_), label = after_stat(p.signif))) +
        theme(
          axis.text.x = element_text(angle = x_angle, hjust = 1, vjust = 1, size = 10),
          legend.position = "bottom",
          legend.title = element_blank()
        )
    }
  )
  return(p4)
}

#' @title 函数：根据值查找变量名
#' @description 比如在函数f(x = a,...)中，运行find_variable_name(x),返回"a"。
#' @param value 输入数据，
#' @return 字符串变量
#' @author *WYK*
find_variable_name <- function(value, env = .GlobalEnv) {
  # 获取环境中的所有变量名
  var_names <- ls(env)

  # 遍历变量名，查找与给定值匹配的变量
  for (var_name in var_names) {
    if (identical(get(var_name, envir = env), value)) {
      return(var_name)
    }
  }

  return(NULL) # 如果没有找到匹配的变量名，返回NULL
}


#'
#' @title 专用Segment plot
#' @param input 绘图数据框
#' @param x_ x轴所在列
#' @param y_ y轴所在列
#' @param color_ 点大小所在列
#' @return ggplot对象
#'
SegmentPoint_p <- function(input, x_, color_, y_) {
  library(tidyverse)

  p1 <- ggplot(data = input, aes(get0(x_), reorder(get0(y_), get0(x_)))) +
    geom_segment(aes_string(xend = 0, yend = y_), linetype = 2) +
    geom_point(aes_string(color = color_), size = 2.75) +
    #   scale_color_continuous(range =c(2,8)) +
    #   scale_x_reverse(breaks = c(0, -0.3, -0.5),
    #                   expand = expansion(mult = c(0.01,.1))) + #左右留空
    ggpubr::theme_pubr() +
    labs(x = get0(color_), y = "", size = "-log10(P)") +
    theme(
      legend.position = "bottom",
      axis.ticks.y.right = element_blank(),
      plot.margin = unit(c(0.2, 0.2, 0.2, 0), "cm")
    )

  return(p1)
}


library(R6)
#'
#' @title 机器学习简单封装的R6类
#' @description ml$new()创建一个新的ml对象，需要三个参数：表达谱df、临床信息df、结果输出路径字符串；
#' - 临床信息df，要求最少有两列，第一列为sample样本信息列，第二列为group，用于样本分组
#' - 举例：ml_obj <- ml$new(
#'    exprs = train_dat$data_exprs[xgene, ],
#'    pheno = train_dat$group,
#'    od = file.path(run_home,'7','2.rf'))
#'    
#' - 设置方法run来进行分析，run中核心参数为method，三选一：method = 'lasso'、'rf'、'svm_rfe'。
#' method = 'lasso'对应lasso分析，非时间依赖性根据输入表达谱以及分组信息筛选重要特征，可选其他参数为seed，设置随机种子，图形自动计算宽高,ml_obj$run(method = "lasso",seed = 1234);
#' method = 'rf'对应随机森林分析，可设置参数：ntree, var_name，seed 三个参数，第一个对应随机森林分析中的ntree数目，默认1000，第二个参数为生成结果的标识字段，默认为'randomforest'，seed 对应分析中的钟子号，默认1110；ml_obj$run(method = "rf",ntree = 1500,var_name = "my_rf")
#' - 可以使用`ml_obj$help()` 来查看简要帮助文档
#' @return 不同方法返回内容不容，所有重要内容均在结果输出路径保存
#'
ml <- R6::R6Class(
  classname = "ml",
  public = list(
    exprs = NULL,
    pheno = NULL,
    od = NULL,
    log = NULL,
    initialize = \(exprs = NULL, pheno = NULL, od = NULL, log = NULL){
      self$exprs <- exprs
      self$pheno <- pheno
      self$od <- od

      if (is.null(log)) {
        self$log <- logging$new("ml")
      }

      index <- list(self$exprs, self$pheno, self$od)

      if (any(map_lgl(index, is.null))) {
        features <- c("expr", "pheno", "od")[which(map_lgl(index, is.null))]
        self$log$warning("缺少{features}信息，需要在`new()`中重新定义")
      }
    },
    help = \(...){
      message("
      help文档：
        - 随机森林中，包含ntree, var_name，seed参数，第一个对应随机森林分析中的ntree数目，默认1000，第二个参数为生成结果的标识字段，默认为'randomforest'，seed 对应分析中的钟子号；
        - lasso中包含seed参数，默认1110;
        - svm-rfe默认使用五倍交叉验证；
        - 使用示例：my_ml_obj$run('lasso',seed = 123.1234)
      ")
    },
    run = function(method = NULL, ...) {
      if (is.null(method)) {
        self$log$error("请选择正确的参数，程序暂时不运行任何分析，参考{.arg method = 'lasso|rf|svm_rfe}'")
        self$help()
        return()
      }

      y <- match.arg(method, c("lasso", "rf", "svm_rfe"))

      f <- function(y) {
        ml_func <- switch(y,
          lasso = \(exprs = self$exprs, od = self$od, pheno = self$pheno, seed = 1110, ...){
            # lasso 建模
            library(glmnet)
            set.seed(seed)
            if (any(is.na(exprs[[2]]))) {
              self$log$warning("LASSO: 表达谱有NA，自动去除...")
              exprs <- na.omit(exprs)
            }

            x_tmp <- exprs %>%
              t() %>%
              as.data.frame() %>%
              rownames_to_column("sample")
            dat <- pheno %>% inner_join(x_tmp)

            x <- dat %>%
              dplyr::select(-any_of(colnames(pheno))) %>%
              as.matrix()
            y <- dat %>%
              dplyr::select(any_of(colnames(pheno))) %>%
              select(-sample) %>%
              pull(1)
            fit1 <<- cv.glmnet(x, y, family = "binomial", type.measure = "deviance", alpha = 1)
            # 绘制图片1
            dir.create(file.path(od), recursive = T, showWarnings = F)
            pdf(file = file.path(od, "lasso_cv.glmnet.fit.pdf"), width = 5, height = 4)
            plot(fit1)
            dev.off()
            # 查看模型各变量系数
            fit1_coef_lambda.min <- coef(fit1, s = fit1$lambda.min)
            # 查看系数非0变量
            fit1.min <- fit1_coef_lambda.min[which(fit1_coef_lambda.min != 0), ]
            # 矩阵转换，使上面结果更方便查看
            fit1.min_2 <- matrix(fit1.min, length(fit1.min), 1)
            rownames(fit1.min_2) <- names(fit1.min)
            colnames(fit1.min_2) <- c("coef")
            lasso_res <- fit1.min_2 %>%
              as.data.frame() %>%
              rownames_to_column("gene")
            lasso_res <- lasso_res[-1, ]
            # (Intercept) -98.7319020

            # 绘制图片2
            fit2 <<- glmnet(x, y, family = "binomial")

            pdf(file = file.path(od, "lasso_glmnet.fit.pdf"), width = 5, height = 4)
            plot(fit2, xvar = "lambda", label = TRUE)
            abline(v = log(fit1$lambda.min), col = "black", lty = 2)
            dev.off()

            lasso_fits <- list(x = x, y = y, fit = fit2, cv.fit = fit1)
            saveRDS(lasso_fits, file.path(od, "lasso_fits.rds"))

            # 输出结果
            write_tsv(lasso_res, file.path(od, "lasso_res.txt"))

            self$log$info("LASSO binomial 回归分析结束")
            lasso_fits$lasso_res <- lasso_res
            lasso_fits
          },
          rf = function(data = self$exprs, group = self$pheno, od = self$od, 
          ntree = 1000, var_name = "random_forest",seed = 1110, ...) {
            self$log$info("随机森林分析开始")
            data <- na.omit(data)

            data <- data %>%
              t() %>%
              as.data.frame() %>%
              rownames_to_column("sample") %>%
              inner_join(group)

            group_col <- colnames(group)[2]

            data[[group_col]] <- factor(data[[group_col]])
            data <- data[, -1]

            dir.create(od, recursive = T, showWarnings = F)

            colnames(data) <- colnames(data) %>%
              str_replace_all("-", "_")

            library(randomForest)
            library(caret)
            set.seed(seed)

            mtry_range <- 1:(ncol(data) / 2)

            # 定义训练控制参数，使用 5 折交叉验证
            train_control <- trainControl(method = "cv", number = 5)

            # 定义 mtry 的候选值
            tune_grid <- expand.grid(.mtry = mtry_range)

            self$log$info("5折交叉验证进行mtry超参数调优")
            rf_model <- train(
              y = data[[group_col]],
              x = data %>%
                select(-any_of(group_col)),
              data = data,
              method = "rf",
              trControl = train_control,
              tuneGrid = tune_grid, importance = T, ntree = ntree
            )

            # bestTune = rf_model$bestTune
            # self$log$info("mtry最优数值：{bestTune}")

            a <- ggplot(rf_model)
            a$data$type <- "normal"
            a$data$type[which.max(a$data[[1]])[1]] <- "best"

            p <- a + geom_point(aes(color = type), show.legend = F) +
              scale_color_manual(values = c("best" = "#ca3434")) +
              egg::theme_article()

            pdf(paste0(od, "/01.", var_name, "_mtry准确率率波动曲线.pdf"), 4.5, 4)
            print(p)
            dev.off()

            pdf(paste0(od, "/02.", var_name, "_gini_准确度.pdf"), 6, 7.4)
            varImpPlot(rf_model$finalModel, main = NULL)
            dev.off()


            pdf(paste0(od, "/03.", var_name, "_错误率波动曲线.pdf"), 5.6, 4.5)
            layout(matrix(c(1, 2), nrow = 1),
              width = c(4, 1)
            )
            par(mar = c(5, 4, 4, 0)) # No margin on the right side
            plot(rf_model$finalModel, main = NULL,lwd= 2,lty=1)
            par(mar = c(5, 0, 4, 2)) # No margin on the left side
            plot(c(0, 1), type = "n", axes = F, xlab = "", ylab = "")
            legend("top", colnames(rf_model$finalModel$err.rate), col = 1:4, cex = 0.8, fill = 1:4)
            dev.off()

            # rf <- randomForest(
            #     y = data[[group_col]], x = data %>%
            #         select(-any_of(group_col)), importance = TRUE, proximity = TRUE, ntree = ntree,
            #     mtry = rf_model$bestTune
            # )
            d <- as.data.frame(importance(rf_model$finalModel)) %>% arrange(desc(MeanDecreaseGini))
            saveRDS(rf_model, file = paste0(od, "/04.", var_name, "_随机森林分析结果.rds"))
            saveRDS(d, file = paste0(od, "/05.", var_name, "_随机森林分析结果_importance_df.rds"))


            self$log$info("随机森林分析结束")
            return(list("d" = d, "train_rf_obj" = rf_model, "best_mtry" = rf_model$bestTune))
          },
          svm_rfe = \(exprs = self$exprs, od = self$od, pheno = self$pheno, seed = 111213,...){
            # SVM
            set.seed(seed)
            features <- exprs %>%
              t() %>%
              as_tibble(rownames = "sample")

            labels <- pheno %>% as_tibble()

            exp_t <- inner_join(features, labels, by = "sample") %>%
              column_to_rownames("sample")

            exp_t[is.na(exp_t)] <- 0

            x <- as.matrix(exp_t[, 1:(ncol(exp_t) - 1)])
            # x <- exp_t
            y <- factor(exp_t[, ncol(exp_t)])

            library(caret)
            library(kernlab)

            self$log$info("SVM-RFE 10x cross validate...")

            control <- rfeControl(functions = caretFuncs, method = "cv", number = 5)
            results <<- rfe(as.matrix(x),
              y = y,
              sizes = seq(1, ncol(x)/2),
              rfeControl = control,
              method = "svmRadial",
              allowParallel = T
            )
            ### 输出结果
            svm_res <- predictors(results)

            self$log$info("SVM-RFE 结束")

            dir.create(od,recursive = F,showWarnings = F)

            write_tsv(tibble(svm_res), file = file.path(od, "svm_res.txt"))
            # 绘制结果
            pdf(file.path(od, "svm_res.pdf"), height = 4, width = 4)
            plot(results, type = c("g", "o"))
            dev.off()

            saveRDS(results,file.path(od, "svm_results.rds"))

            return(results)
          }
        )
        return(ml_func)
      }

      do.call(f(y), list(...))
      # tryCatch(f(),error = \(e){
      #   message('运行失败，请检查。')
      # })
    }
  )
)

#'
#' @title 日志类，生成一个日志对象
#' @description 运用在ml内部显示或生成日志，非重要部分,同样可以通过对象调用`help()`方法查看帮助文档
#' @return NULL
#'
logging <- R6::R6Class("logging",
  public = list(
    name = NULL,
    log_file_path = NULL,
    set_level = NULL,
    set_time_format = NULL,
    set_msg_format = NULL,
    msgs = list(),
    initialize = \(name = "my_log", log_file_path = NULL, set_level = "INFO",
      set_time_format = "%Y-%m-%d %H:%M:%S", set_msg_format = "{level} {name} [{time}]: {msg}"){
      library(stringr)

      self$name <- name
      self$log_file_path <- log_file_path
      if (!is.null(log_file_path)) self$log_file_path <- log_file_path
      self$set_level <- set_level
      self$set_time_format <- set_time_format
      self$set_msg_format <- set_msg_format
    },
    info = \(msg = NULL) {
      msg_ <- private$msg2chr(msg)
      cli::cli_alert_info(msg_)
      self$msgs <- private$collect_msgs(msg_)
      private$save_log()
    },
    warning = \(msg = NULL) {
      msg_ <- private$msg2chr(msg)
      cli::cli_alert_warning(msg_)
      self$msgs <- private$collect_msgs(msg_)
      private$save_log()
    },
    error = \(msg = NULL) {
      msg_ <- private$msg2chr(msg)
      cli::cli_alert_danger(msg_)
      self$msgs <- private$collect_msgs(msg_)
      private$save_log()
    },
    help = \(){
      txt = "
        Create a new logging object using the logging$new() function, with optional parameters:
            - `name` :indicates the log name in log. The default value is 'my_log'
            - `log_file_path`: specifies the log file output path. The default value is `NULL`. No file output is performed
            - `set_level`: indicates the log level. The default value is 'INFO'
            - `set_time_format`:indicates the time format. The default value is '%Y-%m-%d %H:%M:%S'.
            - `set_msg_format`: log format, defaults to '{level} {name} [{time}]: {msg}'

        Usage:
            logger = logging$new(name = 'test',log_file_path = './log_test.log')
            logger$info('test info')
            # ℹ INFO-test [2024-07-19 13:51:58]: test info
            logger$warning('test warning')
            # ! INFO test [2024-07-19 13:51:58]: test warning
            logger$error('test error')
            # ✖ INFO test [2024-07-19 13:51:58]: test error

            # file log_test.log:
            # INFO test [2024-07-19 13:51:58]: test info
            # INFO test [2024-07-19 13:51:58]: test warning
            # INFO test [2024-07-19 13:51:58]: test error

        Other:
            logger$msgs : print all logs
        "
      message(txt)
    }
  ), private = list(
    msg2chr = \(msg_in_msg2chr = NULL){
      formatted_time <- format(Sys.time(), self$set_time_format)

      level <- self$set_level
      name <- self$name
      time <- formatted_time
      msg <- str_glue(msg_in_msg2chr)

      msg_chr <- str_glue(self$set_msg_format)
      return(msg_chr)
    },
    collect_msgs = \(...){
      self$msgs <- c(self$msgs, list(...))
    },
    save_log = \(){
      if (!is.null(self$log_file_path)) {
        if (!file.exists(self$log_file_path)) {
          dir.create(dirname(self$log_file_path), recursive = T, showWarnings = F)
        }
        msg_last <- self$msgs[[length(self$msgs)]]
        readr::write_lines(x = msg_last, file = self$log_file_path, append = TRUE)
      }
    }
  )
)

#'
#' @title 简要nmf分析
#' @param dat 用于nmf分析的表达谱，一般为特征在行，样本在列的df，根据所有特征对样本进行分类
#' @param od 结果输出路径，中间所有结果保存路径
#' @param ranks nmf尝试聚类的聚类数
#' @param best_rank 最终聚类数目，可以先设置一个2或者3然后根据前面跑出来的结果来选一个合适，也可以强行指定聚类数目为2，就固定的聚为2类，但是热图可能没那么好看
#' @param seed 程序运行种子
#' 
#' @return NULL
#'
nmf_simple <- \(dat = tcga_dfs$tumor_exprs[xgene, ], od = "out/1.major_dfs/nmf/",
    ranks = 2:10, best_rank = 3, seed = 1110){

    library(NMF)
    type = 'raw'

    nmf_exprs <- as.matrix(dat)
    stopifnot(min(nmf_exprs) >= 0)

    nmf_res_for_cls <- NMF::nmf(nmf_exprs, ranks, nrun = 50, method = "brunet", seed = seed)
    dir.create(od, recursive = T, showWarnings = F)
    saveRDS(nmf_res_for_cls, file.path(od, "0.nmf根据ranks尝试聚类.rds"))

    nmf_res_for_cls <<- readRDS(file.path(od, "0.nmf根据ranks尝试聚类.rds"))
    pdf(file = file.path(od, "1.0.nmf_progress.pdf"), width = 7.5, height = 6.7,onefile = F)
    show(plot(nmf_res_for_cls))
    dev.off()

    pdf(file = file.path(od, "1.1.nmf_progress_sample_heatmap.pdf"), width = 10, height = 10,onefile = F)
    consensusmap(nmf_res_for_cls,cexRow = 0,cexCol = 0)
    dev.off()

    # 判断最佳rank值的准则就是，cophenetic 值随K变化的最大变动的前点，如3-4变化最大，所以选择最佳rank值为3.
    coph <- nmf_res_for_cls$measures$cophenetic
    coph_diff <- NULL
    for (i in 2:length(coph)) {
        coph_diff <- c(coph_diff, coph[i - 1] - coph[i])
    }
    k.best <- which.max(coph_diff) + 1

    cli::cli_alert_info(cli::style_bold(">> 自动计算best_rank {k.best}，程序运行设置为 {best_rank}"))

    # 再次NMF,rank=3
    nmf_best_rank <- nmf(nmf_exprs,
        rank = best_rank,
        nrun = 50,
        seed = seed,
        method = "brunet"
    )

    saveRDS(nmf_best_rank, file.path(od, "2.0.选用最佳聚类数后的nmf聚类结果.rds"))

    # 常见的 NMF 方法
    # Brunet (brunet):    
    # 描述: 使用多重更新规则（Multiplicative Update Rules）。
    # 应用场景: 适用于一般的 NMF 分解，尤其是基因表达数据的聚类分析。
    # 优点: 简单且易于实现，收敛速度较快。

    # Lee & Seung (lee):    
    # 描述: 采用 Lee 和 Seung 提出的乘法更新规则。
    # 应用场景: 适用于图像处理和文本挖掘等领域。
    # 优点: 理论基础扎实，广泛应用于各种数据类型。

    # KL Divergence (kl):    
    # 描述: 基于 Kullback-Leibler 散度的优化方法。
    # 应用场景: 适用于概率分布数据的分解，如文本数据的主题模型。
    # 优点: 对稀疏数据表现良好，适合处理概率分布。

    # Frobenius Norm (frobenius):    
    # 描述: 最小化 Frobenius 范数的优化方法。
    # 应用场景: 适用于一般的矩阵分解问题。
    # 优点: 简单且计算效率高。

    # Offset (offset):    
    # 描述: 在 NMF 分解中加入偏移量的优化方法。
    # 应用场景: 适用于数据中存在偏移量的情况，如基因表达数据。
    # 优点: 能够处理数据中的偏移量，提高分解的准确性。

    res_methods <- nmf(nmf_exprs, best_rank, list("lee", "brunet", "nsNMF"), nrun = 50, seed = seed)
    pdf(
      file = file.path(od, "2.1.选用最佳聚类数后的nmf聚类，方法选择heatmap，默认brunet.pdf"),
      width = 6.5, height = 6.5,onefile = F
    )
    consensusmap(res_methods, cexRow = 0, cexCol = 0)
    dev.off()

    # 可视化结果，样本聚类热图
    pdf(file.path(od, "2.2.选用最佳聚类数以及brunet后，未提取特征的nmd聚类fig.pdf"), width = 5, height = 4.5,onefile = F)
    consensusmap(nmf_best_rank,
        labRow = NA,
        labCol = NA
    )
    dev.off()

    # 提取特征
    index <- extractFeatures(nmf_best_rank, "max")
    sig.order <- unlist(index) %>% na.omit()
    nmf.input2 <- nmf_exprs[sig.order, ]

    rownames(nmf_exprs)[sig.order]

    saveRDS(rownames(nmf_exprs)[sig.order], file.path(od, "3.重要特征.rds"))

    # 用特征再做NMF
    nmf_best_rank_res <- nmf(nmf.input2,
        rank = best_rank,
        seed = seed,
        method = "brunet"
    )
    saveRDS(nmf_best_rank_res, file.path(od, "4.选用最佳聚类数以及提取重要特征后的nmf聚类结果.rds"))


    # 设置颜色
    jco <- ggsci::pal_d3()(best_rank)

    switch(type,
      raw = {
        group <- predict(nmf_best_rank) # 提出亚型

        d <- as.data.frame(group) %>%
          rownames_to_column("sample") %>%
          mutate(group = str_c("C", group)) %>%
          rename(nmf_cls = 2)
        write_tsv(d, file.path(od, "5.0.nmf_sample_cls_raw_features.tsv"))
      },
      selected = {
        group <- predict(nmf_best_rank_res) # 提出亚型

        d <- as.data.frame(group) %>%
          rownames_to_column("sample") %>%
          mutate(group = str_c("C", group)) %>%
          rename(nmf_cls = 2)
        write_tsv(d, file.path(od, "5.1.nmf_sample_cls_with_nmf_selected_feature.tsv"))
      }
    )


    

    # 可视化结果，样本聚类热图
    pdf(file.path(od, "6.选用最佳聚类数以及提取重要特征后的样本聚类heatmap.pdf"), width = 5, height = 4.5,onefile = F)
    consensusmap(nmf_best_rank_res,
        labRow = NA,
        labCol = NA,
        annCol = data.frame("cluster" = group[colnames(nmf_exprs)]),
        annColors = list(cluster = c("1" = jco[1], "2" = jco[2], "3" = jco[3]))
    )
    dev.off()

    names(jco) <- seq_len(best_rank)
    pdf(file.path(od, "7.选用最佳聚类数后特征heatmap.pdf"), width = 5, height = 4.5,onefile = F)
    basismap(nmf_best_rank,
        cexCol = 1,
        cexRow = 1,
        annColors = list(jco)
    )
    dev.off()
}


#'
#' @title 检查文件路径的文件夹是否存在
#' @param x 文件路径
#' @return x
#'
file_dir_check <- check_file_dir <- \(x = NULL){
  dir.create(dirname(str_glue(x)), showWarnings = F, recursive = T)
  return(str_glue(x))
}


#'
#' @title 从网页HaploRegv4.2 在线判断该snp是否具备eQTL信息
#' @param snp 字符串，snp id，eg：rs1800440
#' @return 逻辑值，TURE or FALSE
#'
get_snp_eQTL_in_HaploRegv4.2 <- \(snp){
    library(rvest)
    library(httr)
    url <- "https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php"

    logger = logging$new(name = "get_snp_eQTL_in_HaploRegv4.2",set_msg_format = "{name} [{time}]: {msg}")
    logger$info(snp)    
    # 创建POST请求的表单数据
    form_data <- list(
        query = snp,
        gwas = "",
        ldThresh = "0.8",
        ldPop = "EUR",
        epi = "vanilla",
        cons = "siphy",
        genetypes = "gencode",
        snpfunc = "vanilla",
        dbSNP = "150",
        submit = "submit"
    )

    # 发送POST请求
    response <- POST(url, body = form_data, encode = "form")

    # 检查请求是否成功
    if (status_code(response) == 200) {
        # 获取原始内容并指定编码
        raw_content <- httr::content(response, as = "raw")
        page <- read_html(raw_content, encoding = "ISO-8859-1")

        # 提取感兴趣的内容（例如表格数据）
        results <- page %>%
            html_nodes("table") %>%
            html_table(fill = TRUE)
    } else {
        logger$warning("请求失败，状态码：{status_code(response)}, snp: {snp}")
    }

    if(length(results) < 4){
        a = FALSE
        names(a) = snp
        return(FALSE)
    }

    res_table <- results[[4]] %>% as.data.frame()

    colnames(res_table) <- res_table[1, ] %>% as.character()
    res_table <- res_table[-1, ]

    eQTLs <- res_table %>%
        filter(variant == snp)

    if (eQTLs[["Selected eQTLhits"]] == "") {
        a = FALSE
        names(a) = snp
        return(a)
    } else {
        a = TRUE
        names(a) = snp
        return(a)
    }
}


#'
#' @title 一对多相关性分析，分别使用dat1的特征与dat2的所有特征进行相关性分析
#' @param dat1 特征在行，样本在列的dataframe
#' @param dat2 特征在行，样本在列的dataframe
#' @param od 结果输出路径
#' @param w 图宽
#' @param h 图高
#' @return NULL
one_with_more_cor <- \(dat1, dat2, method = "spearman", od, w = 3, h = 4){

  dir.create(od,recursive = T,showWarnings = F)
  common_smaples <- intersect(colnames(dat1), colnames(dat2))
  common_g <- rownames(dat1)

  dat1_ready <- dat1[, common_smaples]
  dat2_ready <- dat2[, common_smaples]

  hub_immune_cor <- WGCNA::corAndPvalue(t(dat1_ready) %>% as.matrix(),
    dat2_ready %>% t() %>% as.matrix(),
    method = "spearman"
  )
  hub_immune_cor_dat <- data.frame(hub_immune_cor$cor, check.names = F) %>%
    rownames_to_column(var = "xname") %>%
    pivot_longer(cols = -xname, values_to = "cor", names_to = "yname")

  hub_immune_dat <- data.frame(hub_immune_cor$p, check.names = F) %>%
    rownames_to_column(var = "xname") %>%
    pivot_longer(cols = -xname, values_to = "pvalue", names_to = "yname") %>%
    mutate(cor = hub_immune_cor_dat$cor, "-log10pvalue" = -log10(pvalue))

  write.table(hub_immune_dat, str_glue("{od}/dat1_dat2_correlation.txt"), row.names = F, sep = "\t", quote = F)

  plot_list <- lapply(common_g, function(x) {
    x_dat <- hub_immune_dat %>%
      filter(xname == x) %>%
      arrange(desc(cor))
      
    x_dat$ynum <- 1:nrow(x_dat)

    x_plot <- ggplot(x_dat, aes(x = ynum, y = cor)) +
      geom_segment(aes(x = ynum, xend = ynum, y = 0, yend = cor), color = "gray") +
      geom_point(aes(color = pvalue, size = cor)) +
      scale_color_viridis_c(direction = -1) +
      scale_x_continuous(
        breaks = x_dat$ynum, labels = x_dat$yname,
        sec.axis = sec_axis(~., breaks = c(x_dat$ynum, nrow(x_dat) + 1), labels = c(round(x_dat$pvalue, 3), "pvalue"))
      ) +
      scale_size(range = c(2, 5)) +
      labs(y = "Correlation", x = NULL, title = x) +
      coord_flip() +
      theme_bw() +
      theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5), axis.ticks.y.right = element_blank())
    ggsave(str_glue("{od}/{x}_cor.pdf"), x_plot, width = w, height = h)
    return(x_plot)
  })
  return(plot_list)
}




#'
#' @title 利用ncbi网页爬虫目标snp的Alleles类型
#' @param rs_id 输入rs_id
#' @return 字符串，该位点可能突变的碱基情况
get_rs_alleles <- \(rs_id){
    # 设置目标URL
    url <- str_glue("https://www.ncbi.nlm.nih.gov/snp/{rs_id}")

    library(rvest)
    # 读取网页内容
    webpage <- read_html(url)

    a <- webpage %>%
        html_nodes("body") %>%
        html_text()
    a %<>% str_split("\n") %>% unlist()

    a <- a[which(grepl("Alleles", a)) + 4] %>% gsub("[^A-Za-z]", "", .)

    a %<>%
        str_split_1("") %>%
        unique() %>%
        str_c(collapse = "/")

    return(str_glue("[{a}]"))
}

#'
#' @title 利用ncbi网页爬虫目标snp的前后5000碱基序列
#' @param rs_id 输入rs_id
#' @param gap 前后序列长度，默认5000
#' @return fasta格式的字符串
get_snp_fasta <- \(rs_id, gap = 5000){
    rs_id_input <- rs_id
    # logger <- logging$new('get_snp_fasta',set_msg_format = '{name} [{time}]: {msg}')

    if (!exists("rmvar", envir = .GlobalEnv)) {
        rmvar <<- data.table::fread("/Pub/Users/fuxj/Database/RMVar/RMVar_Human_basic_info_m6A.txt", sep = "\t", header = T, fill = T)
    }

    rs_df <- rmvar %>%
        filter(rs_id %in% rs_id_input) %>%
        select(chromosome, position, gene, rs_id, strand) %>%
        distinct()

    loop_index <- seq_len(nrow(rs_df))
    get_snp_fasta_inner <- \(x){
        rs_df <- as.data.frame(rs_df)
        chr <- rs_df[x, "chromosome"]
        pos <- rs_df[x, "position"]
        # logger$info("get {rs_df[x, 'rs_id']}'s allele")
        # allele <- get_rs_alleles(rs_df[x, "rs_id"]) %>% str_split_1('') %>% .[2]
        offset <- gap
        type <- rs_df[x, "strand"]

        library(BSgenome)
        library(BSgenome.Hsapiens.UCSC.hg19)

        seq <- paste0(
            getSeq(Hsapiens, chr, start = (pos - offset), end = (pos - 1), strand = type),
            getSeq(Hsapiens, chr, start = pos, end = pos, strand = type),
            getSeq(Hsapiens, chr, start = (pos + 1), end = (pos + offset), strand = type)
        )
        return(str_glue('>{rs_df[x, "rs_id"]}\n{seq}\n'))
    }
    res <- map(loop_index, get_snp_fasta_inner)
    res[[1]]
}


#'
#' @title 利用ncbi网页爬虫目标snp的前后5000碱基序列
#' @param symbol 基因SYMBOL
#' @return data.frame 
get_gene_entriz_df <- function(symbol) {
    library(org.Hs.eg.db)
    # 使用rentrez的esearch功能找到基因ID
    if (!exists(".gene_id_db",envir = .GlobalEnv)) {
        .gene_id_db <<- AnnotationDbi::select(org.Hs.eg.db,
            keys = keys(org.Hs.eg.db, keytype = "SYMBOL"),
            columns = c("SYMBOL", "ENTREZID"), keytype = "SYMBOL"
        )
    }

    search_result <- .gene_id_db %>%
        filter(SYMBOL == symbol) 

    search_result[1,]
}


#'
#' @title 利用种植处理重复表达谱
#' @param expr 表达谱
#' @param feature_col 特征所在列，比如说基因symbol
#' @param type 重复特征标准化类型，默认为median
#' @return data.frame 
deal_dup_expr <- \(expr, feature_col = "gene", type = "median"){
  colnames(expr)[colnames(expr) == feature_col] <- "gene"

  which(duplicated(expr$gene)) -> i

  cli::cli_alert_info("{length(i)} features is duplicated in exprs.")

  expr %>% filter(gene %in% expr$gene[i]) -> dup_expr
  expr %>% filter(!gene %in% expr$gene[i]) -> no_dup_expr

  dup_expr %<>% group_by(gene) %>% summarise_all(.funs = get0(type))

  expr <- bind_rows(no_dup_expr, dup_expr) %>%
    as.data.frame() %>%
    na.omit()

  tibble::column_to_rownames(expr, "gene") -> expr
}